#installaation de package nécéssaire en local
#install.packages("BiocManager")
#BiocManager::install("phyloseq")
#install.packages("phyloseq")
#install.packages("remotes")
#remotes::install_github("mahendra-mariadassou/phyloseq-extended", ref = "dev")
#install.packages("gridExtra")
#Installation des librairies
library(tidyverse)
library(phyloseq)
library(phyloseq.extended)
library(ggplot2) ## for graphs necessary ???
library(ape) ## for tree manipulation
library(vegan) ## for community ecology analyses
library(scales)
library(grid)
library(gridExtra) ## for design graph
library(knitr)
library(kableExtra)
chemin <- "~/save/Milint/Resultat/"
load("BIOM.Rdata")Le BIOM a été créé dans le script01. Il est cependant brut
Somme de tous les reads et on applique un filtre de 10-5 afin de s’assurer que ce ne sont pas des artefacts Dans notre cas permet de filtrer les taxa avec des valeurs supérieures à 180 reads par OTU
[1] 5042 240312
[1] 5003 239715
# moyenne
moy <- sum(taxa_sums(Biom))/846
moyf <- sum(taxa_sums(filter_biom))/846
mi <- min(taxa_sums(Biom))
ma <- max(taxa_sums(Biom))
mif <- min(taxa_sums(filter_biom))
maf <- max(taxa_sums(filter_biom))
med <- median(taxa_sums(Biom))
medf <- median(taxa_sums(filter_biom))
df <- data.frame(
échantillon = c(846,846),
mean = c(moy,moyf),
median = c(med,medf),
min = c(mi,mif),
max = c(ma, maf)
)
rownames(df) <- c("Biom","Filter_Biom")
kable(df)| échantillon | mean | median | min | max | |
|---|---|---|---|---|---|
| Biom | 846 | 21339.08 | 5 | 0 | 2709474 |
| Filter_Biom | 846 | 21207.16 | 1092 | 181 | 2709474 |
Le niveau d’abondance semble hétérogène. Il faudra donc prévoir un procédé de rarefaction
plot_bar(Biom)+
ggtitle("Abondance par échantillon sur données totales")+
theme(axis.text.x=element_blank(),
axis.ticks.x=element_blank())plot_bar(filter_biom)+
ggtitle("Abondance par échantillon sur données filtrées")+
theme(axis.text.x=element_blank(),
axis.ticks.x=element_blank())Dorénavant nous ne travaillerons que sur données filtrées Ne sert à rien d’aller au niveau Phylum ou infra car trop d’échantillons donc aucune couleur n’est vissible Nous remarquons que le niveau d’abondance par échantillon est hétérogène de 25 000 en moyenne à 250 000.
J’avais fait des niveaux d’abondance par échantillon en fonction de différentes variables (sexe, satut tabagique…) mais ne sert à rien car ça équivaut à voir la profondeur et ce n’est pas biologiquement informatif donc enlevé.
| Var1 | Freq |
|---|---|
| Actinobacteria | 22 |
| Bacteroidetes | 223 |
| Firmicutes | 686 |
| Fusobacteria | 2 |
| Proteobacteria | 45 |
| Synergistetes | 9 |
| Tenericutes | 5 |
| Verrucomicrobia | 5 |
# Coloration ne fonctionne pas
# scale_fill_manual( values = c("Actinobacteria" = "black", "Bacteroidetes" = "orange", "Firmicutes" = "blue",
# Fusobacteria = "yellow",
# Proteobacteria = "red",
# Synergistetes= "pink",
# Tenericutes = "green",
# Verrucomicrobia = "purple"))
#
#geom_bar(aes((x="Phylum) colour = "Phylum")#Abondance par phylum
phy <- tax_glom(filter_biom, taxrank = "Phylum")
plotdata<-psmelt(phy)
ggplot(plotdata,aes(y=Abundance))+
geom_boxplot(aes(group="Phylum", fill="Phylum"))+
facet_wrap(~Phylum, scales="free_y",nrow=2)plot_composition(filter_biom, "Phylum", "Firmicutes", "Family", fill = "Family")+
facet_grid(scales = "free_x") +
theme(axis.text.x = element_blank())plot_composition(filter_biom, "Phylum", "Firmicutes", "Genus", fill = "Genus", numberOfTaxa = 20) +
facet_grid(scales = "free_x", space = "free_x") +
theme(axis.text.x = element_blank())plot_composition(filter_biom, "Phylum", "Firmicutes", "Genus", fill = "Genus", numberOfTaxa = 5) +
facet_grid(~tabac, scales = "free_x", space = "free_x")+
theme(axis.text.x = element_blank())plot_composition(filter_biom, "Phylum", "Firmicutes", "Genus", fill = "Genus", numberOfTaxa = 5) +
facet_grid(~AGE, scales = "free_x", space = "free_x")+
theme(axis.text.x = element_blank()) * Composition en Firmicutes au sein du groupe sexe
plot_composition(filter_biom, "Phylum", "Firmicutes", "Genus", fill = "Genus", numberOfTaxa = 5) +
facet_grid(~SEX, scales = "free_x", space = "free_x")+
theme(axis.text.x = element_blank())Au vu des différents graphs, il y a une grande diversité et nous ne voyons pas de genre majoritaire quelque soit la variable. Pour l’ordre Clostrdiales est majoritaires quelque soit la variable prise (exemple ci dessous)
plot_composition(filter_biom, "Phylum", "Firmicutes","Order", fill = "Order") +
facet_grid(~AGE, scales = "free_x", space = "free_x")+
theme(axis.text.x = element_blank())Je pense que les graphs suivants ne sont d’aucune utilité . J’ai fait les mêmes sur les femmes mais trop d’échantillons
knitr::opts_chunk$set(
fig.width = 6, fig.height = 15,
fig.align = "center")
# échantillon ayant plus de 50000 abondance
#prune_samples(sample_sums(biom) > 50000, biom)
#Abondance sur Hommes ayant plus de 25 000
H_more_25K <- subset_samples(filter_biom, SEX == 1 & sample_sums(filter_biom) > 25000)
H_m_25K <- plot_bar(H_more_25K, fill = "Phylum")
#Abondance sur Hommes ayant moins de 25 000
H_less_25K <- subset_samples(filter_biom, SEX == 1 & sample_sums(filter_biom) <= 25000)
H_l_25k <- plot_bar(H_less_25K, fill = "Phylum")
#plot_composition(H_less_20K, "Phylum", "Firmicutes", "Family", fill = "Family")
grid.arrange(H_m_25K,H_l_25k , nrow=2,
top=textGrob("abondance sup ou inf à 25K catégorie Hommes ",gp=gpar(fontsize=20,font=3)))Exploration de l’impact de chaque covariable sur la diversité-alpha
plot_richness(filter_rare, x = "AGE", measures = c("Observed", "Shannon", "InvSimpson")) +
geom_boxplot(aes(group = AGE)) ## add one boxplot per AGE valuesAge diffère significativement en terme d’observation d’OTU
div_data <- cbind(estimate_richness(filter_biom), ## diversity indices
sample_data(filter_biom) ## covariates
)
model <- aov(Observed ~ AGE, data = div_data)
anova(model)Mais également en terme de diversité en nombre d’espèce
plot_richness(filter_rare, x = "SEX", measures = c("Observed", "Shannon", "InvSimpson")) +
geom_boxplot(aes(group = SEX)) ## add one boxplot per SEX valuesPas de différence signicative pour le nbre d’OTU mais différence significative pour le nbre d’espèces
plot_richness(filter_rare, x = "tabac", measures = c("Observed", "Shannon", "InvSimpson")) +
geom_boxplot(aes(group = tabac)) ## add one boxplot per AGE values